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The origin of the vortex-core states in s-wave and d x 2_ y 2-~w&ve superconductors is investigated 
by means of some selected numerical experiments. By relaxing the self-consistency condition in the 
Bogoliubov-de Gennes equations and tuning the order parameter in the core region, it is shown that 
the suppression of the superfluid density in the core is not a necessary condition for the core states 
to form. This excludes "potential well" types of interpretations for the core states. The topological 
defect in the phase of the order parameter, however, plays a crucial role. This observation is 
explained by considering the effect of the vortex supercurrent on the Bogoliubov quasiparticles, 
and illustrated by comparing conventional vortices with multiply-quantized vortices and vortex- 
antivortex pairs. The core states are also found to be extremely robust against random phase 
disorder. 
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I. INTRODUCTION 

The vortices govern the electromagnetic response of 
type-II superconductors and have been extensively stud- 
ied, both experimentally and theoretically. 1,2 A vortex is 
formed from a core of radius r c ~ £ where the superfluid 
density is gradually suppressed, £ the superconducting 
coherence length, and it is surrounded by a supercurrent 
which screens the magnetic field on a length of order 
A, the penetration depth. The vortices are strong inho- 
mogeneities of the superconducting condensate and they 
scatter the quasiparticles in several different ways. 3 In 
particular, the vortices can capture Bogoliubov excita- 
tions into low-energy localized states. 

The vortex-core states play an important role in the 
thermodynamic and transport properties in the mixed 
state. For example, when vortices move in an applied 
electric field, the core states interact with the lattice 
and are thus responsible for the dissipation of energy. 
These states are also affected by localized perturbations 
and they contribute to the pinning of vortices by de- 
fects or impurities. Recently, mesoscopic superconduct- 
ing disks have attracted much attention. 4 In these sys- 
tems a "giant-vortex state" can be stabilized, where a 
single vortex at the center of the sample carries the whole 
magnetic flux. In such a case the core states are the main 
low-energy excitations, and they are expected to play a 
dominant role. Furthermore, the strong dependence of 
the vortex-core energy spectrum on the applied magnetic 
field opens interesting perspectives for applications. 5 

In s-wave superconductors the vortex-core bound 
states were predicted long ago, based on approximate so- 
lutions of the microscopic Bogoliubov-de Gennes (BdG) 
equations, 6,7 and subsequently observed in NbSe2 using 
scanning tunneling spectroscopy. 8 The early analytical 
results were confirmed by a complete numerical solution 
of the BdG equations. 9 Extended quasiparticle excita- 
tions in the mixed state were often studied within the 
quasiclassical approach. 2,10 Although this approximation 
is considered inaccurate near the vortex core, it was used 
by Volovik to argue that the number of branches of core 



states crossing the Fermi energy as a function of angular 
momentum (considered as a continuous variable) is equal 
to the winding number of the vortex. 11 This prediction 
was confirmed recently by a detailed numerical solution 
of the BdG equations for multiply-quantized vortices. 12 
The vortex-core states are usually thought of as Andreev 
bound states, i.e. standing waves resulting from the mul- 
tiple Andreev reflection at the normal/superconductor 
boundary in the core. 2,10,13,14 This interpretation sug- 
gests that the suppression of the superfluid density in 
the core is the main reason for the formation of the 
core states. The vanishing of the superfluid density in 
the core, however, plays no role in Volovik's argument. 
Instead, the structure of the vortex-core energy spec- 
trum in the approach of Ref. 11 is entirely determined 
by the winding number of the vortex, which measures 
the strength of the supercurrent circulating around it. 
This result seems difficult to reconcile with the Andreev- 
bound state picture. In particular, the peculiar depen- 
dence of the spectrum of core states on vorticity 12 can 
hardly be attributed to the order-parameter suppression 
alone. Hence it is of interest to identify the roles played 
by the supercurrent, on the one hand, and by the order- 
parameter suppression, on the other hand, in the forma- 
tion of the vortex-core states. 



Based on numerical and analytical solutions of the mi- 
croscopic BdG equations, I show that the structure of 
the bound-state spectrum in s-wave and c£-wave vortices 
is determined by a topological constraint that the circu- 
lating superfluid imposes to Bogoliubov quasiparticles. 
The suppression of the order parameter in the core plays 
a minor quantitative role, slightly changing the energy 
of the states with small angular momenta. This implies, 
in particular, that a complete self-consistent treatment 
of the order-parameter profile is in general not necessary, 
unless one is interested in detailed quantitative predic- 
tions. These results suggest that the mechanism leading 
to quasiparticle localization in vortices is quite different 
from other localization mechanisms in condensed matter. 
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II. TOPOLOGY OF VORTICES AND BDG 
EQUATIONS 

A vortex is a topologically stable defect of the super- 
conducting order parameter ^(r) = A(r)e lx ^ r h 15 It is 
characterized by a winding number u, a topological in- 
variant defined as 
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The integral runs along a closed path around the vortex 
axis. \v\ counts the number of 2-7T rotations of the phase 
x(r) along the path. Clearly v is invariant under all con- 
tinuous distortions of the phase. Except near surfaces or 
interfaces, v corresponds to the magnetic flux $ carried 
by the vortex according to $ = i/$ with $ = ^ the 
superconducting flux quantum. Another characteristics 
of the vortex is the shape of the order-parameter mod- 
ulus A(r) in the vicinity of the vortex axis, where it is 
constrained to vanish. 

For an isolated axisymmetric vortex in a s-wave super- 
conductor one has A(r, 8, z) = A(r) and %(r, 8, z) = v6, 
where v is the winding number consistently with Eq. (1). 
The modulus A(r) vanishes as r — > 0, and approaches 
the constant value Aoo at r > r c , where r c <~ £ <C A 
in type-II superconductors. The excitation spectrum of 
the vortex is determined by the Bogoliubov-de Gennes 
(BdG) equations, 
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where h = ^-(p — eA) 2 — E F and u (v) is the electron 
(hole) wave function of the excitation. In order to solve 
Eq. (2) one usually eliminates the phase x °f the order 
parameter from the off-diagonal terms by performing the 
gauge transformation A — > A — j^^X- The order pa- 
rameter changes according to ^ — > ^e^ lx and thus be- 
comes real. 6 This transformation is achieved by writing 
the wave functions as 
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The phase ±jj 8 in Eq. (3) is an Aharonov-Bohm phase, 
reflecting the fact that the gauge function —\x carries a 
singular magnetic field at the origin. 16 Furthermore, the 
quantum number p must be such that the total phase 
accumulated by the electron and hole upon a 2ir rotation 
around the origin is consistent with the enclosed flux, i.e. 
v 

+ r 
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Inserting Eq. (3) into Eq. (2) solves for the 9 and z de- 
pendencies and leads to the radial equations for the real 
functions ip±: 
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Here g = k F t; and I have introduced the dimensionlcss 
variables p = r/£, e = E/E F , and 5(p) = A(p£)/E F . For 
simplicity I restricted in Eq. (5b) to the two-dimensional 
case by putting k z = 0. The vector potential A was also 
omitted since it is small in the core region compared to 
the gauge field ^Vx when A > £. 6 Indeed A - j^r, 
and the ratio of A to the gauge field is thus of order 
(r/A) 2 . 

Although the gauge transformation removes the phase 
of the order parameter, it does not eliminate the super- 
current from the problem. In the radial equation, the 
supercurrent shows up as a central potential containing a 
repulsive term (v/2p) 2 , as well as a term ±pv/ p 2 which is 
either attractive or repulsive: if v < (supercurrent flow- 
ing counter-clockwise, corresponding to a positive mag- 
netic field along the z axis) this term is attractive (repul- 
sive) for the electrons (holes) which move alike the su- 
pcrfmid (p > 0). Therefore, the supercurrent acts on the 
electron and hole parts of the BdG excitations in differ- 
ent ways, and tends to decrease the angular momentum 
of the former and to increase the angular momentum of 
the latter (for v < and p > 0). This effect of the super- 
current on the BdG excitations is central to understand 
the formation of the core states (see Sect. IIIB). Further- 
more, the strength of the supercurrent fixes the parity of 
the vortex-states angular momentum in Eq. (4), which 
is half an odd integer for odd v and integer for even v. 
When v is odd, the flux carried by the vortex is not a 
multiple of the flux carried by the quasiparticle, hence 
a topological frustration which translates into a branch 
cut discontinuity — removed by the gauge transformation 
in Eq. (3) — in the angular wave function. 

Eq. (2) possess the well-known particle-hole sym- 
metry (u, v, E) <-» (v* , — u* , —E). Looking at the 
wave functions in Eq. (3), one sees that in the ra- 
dial equation this symmetry becomes (p, ip-, s) ^ 
(—p, tp-i —ip+, —£)'■ the vortex-core energy spectrum is 
invariant under the simultaneous inversion of angular mo- 
mentum and energy. Furthermore, it is clear from Eq. (5) 
that the spectrum is also invariant under the simultane- 
ous inversion of p and v. 



III. THE ISOLATED S-WAVE VORTEX 

Self-consistent solutions of the BdG equations for the 
order-parameter profile and the energy spectrum of an 
isolated s-wave vortex were already reported, both in the 
singly-quantized 9 and multiply-quantized 12 cases. The 
purpose of this section is to repeat these calculations 
without achieving the self-consistency in A(r), in order 
to clarify the role played by the detailed form of A(r) 
and by the winding number v in the formation of the 
core states. Analytical solutions for the singly-quantized 
vortex also exist. 6 ' 7 I will discuss a new analytical solu- 
tion which is valid for all (integer) values of v, and which 
emphasizes the key role of the winding number. 
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FIG. 1: Subgap energy spectrum of a singly-quantized s-wave 
vortex as a function of angular momentum u for various order- 
parameter profiles in the core: nearly self-consistent profile 
(black circles), "no-core" profile (white circles), and step- like 
profile (squares). The triangles show the energy spectrum of 
a normal region analogous to a vortex core embedded in a 
uniform superconductor (no supercurrent). 

A. Numerical results 

The BdG equations (5) were solved numerically using 
the Bcsscl-function expansion described in Ref. 9. Be- 
sides its winding number and order-parameter profile, 
an isolated vortex in a continuum free-electron model 
is characterized by the bulk parameters g = /cf£ and 
S = Aoo/Ep. Physical values of g range from <~ 1 in 
high-T c materials to 10-10 4 in conventional supercon- 
ductors. Simulations were performed for g between 1 
and 100 (the computational effort increases rapidly with 
increasing g). The parameters g and 8 are in principle 
related by the BCS relation £ w i.e. gS w 2/tt. For 

each g, values of S between 0.1 and 10 times the BCS 
value were considered. The reported conclusions apply 
to the whole domain of parameters investigated. Below 
I present results for g = 10 and S = 2/(irg). 17 For the 
order-parameter profile the following form was used: 18 

( r \ M 
A(r) = Aoo tanh = . (6) 

\VW\ r cJ 

This functional form with r c = £ is in good qualitative 
agreement with the self-consistent results of Ref. 12. By 
tuning the value of r c one can study the effect of the gap 
profile on the vortex-core spectrum. In particular, the 
limit r c — > corresponds to a vortex with no normal core. 
(In the remainder of this paper, I shall use the expression 
normal core as a synonym for suppression of A(r) in 
the vortex core.) According to the Andreev- reflection 
picture, one may expect that reducing the core radius 
will increase the energy separation between the vortex- 
core energy levels, and eventually push the core states 
outside the energy gap into the continuum as r c — > 0. 

The spectra of subgap states for a singly-quantized 
vortex {v = —1) and for gap profiles corresponding to 
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FIG. 2: Local density of states for the systems shown in Fig. 1: 
conventional singly-quantized vortex (left), vortex without 
normal core (center), and normal core without supercurrent 
(right). A thermal broadening was used with feT = Aoo/10. 



r c = £ and r c — are shown in Fig. 1. The eigen- 
values are displayed as a function of the angular quan- 
tum number /j,. At energies outside the superconduct- 
ing gap, the BdG states form a continuum not shown in 
the figure. The spectrum obtained for a step-like pro- 
file A(r) = Aoo#(r — £), which is often used in analyt- 
ical calculations, is also displayed for comparison. The 
spectra are similar in all three cases, except for small dif- 
ferences at low values of \n\. These differences have a 
rather simple explanation: the corresponding cigenstates 
are concentrated close to the vortex axis — the maximum 
of the wave function lies roughly at position |/i|/.g in units 
of £; thus the pairing energy of the states with < g 
is lowered when the order parameter gets suppressed at 
r < £. It is clear from the figure, however, that the spec- 
trum retains its general shape even when the vortex has 
no normal core. In particular, the number of branches 
crossing the Fermi energy is always one, in agreement 
with Volovik's theorem. 

We have seen that suppressing the normal core, keep- 
ing only the supercurrent, does not change the core states 
qualitatively. We may now do the converse: in order to 
suppress the supercurrent and keep only the core, we set 
the phase of the order parameter to zero [i.e. v = in 
Eqs. (3)-(5)] and we use the profile Eq. (6) with v = 1 
and r c — £. The resulting order parameter no longer 
describes a vortex, since the winding number vanishes, 
but just a normal region embedded in a uniform super- 
conductor. The resulting spectrum (triangles in Fig. 1) 
is qualitatively different from the spectrum of a singly- 
quantized vortex. Consistently with Volovik's result, no 
branch of core states cross the Fermi level, and therefore 
no low-energy states exist, although multiple Andreev 
reflections are in principle possible in this system. The 
suppression of the order parameter slightly decreases the 
pairing energy of the low-|/x| electron and hole excita- 
tions of the uniform superconductor, and gives rise to 
two shallow branches of states near the gap edges. 
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FIG. 3: Energy spectra for multiply-quantized s-wave vor- 
tices. The black circles correspond to vortices having a nor- 
mal core given by Eq. (6), and the white circles correspond 
to vortices without normal core. The dashed vertical lines 
delimit the region \p\ < \u\g where the levels are sensitive to 
the profile of the order parameter. 



Fig. 2 compares the local density of states (LDOS) of 
a singly-quantized vortex with and without normal core. 
As can be inferred from Fig. 1, the effect of the normal 
core is mainly to raise the peak at the vortex center, 
without changing the structure of the LDOS. The energy 
of the peak is also closer to the Fermi energy when the 
normal core is present. This is due to the core-induced 
energy change at low p. Indeed, since only states with 
/u±| = contribute to the LDOS at r = [see Eq. (Al)], 
we have 

JV(0, £0 a £<*(25 -!$=-*) (i/<0), (7) 



where the index i corresponds to the various eigenvalues 
at fixed p. Fig. 2 also shows the LDOS of a normal 
core without supercurrent. In this case there is no zero- 
bias peak, but two sharp peaks near ±Aoo corresponding 
to the two branches in Fig. 1. These peaks disappear 
beyond a distance of order £ and the LDOS becomes that 
of the uniform superconductor for r > £. In contrast, for 
a true vortex the perturbation to the bulk DOS extends 
to distances much larger than £, irrespective of the core 
radius r c , due to the high-|/x| states in Fig. 1. 

The vortex-core spectra and the LDOS of multiply- 
quantized s-wave vortices (\v\ > 1) were calculated in 
Rcf. 12, using the self-consistent order parameter pro- 
file. The number of branches crossing Ey was as 



expected. I have repeated these calculations using the 
profile Eq. (6) with both r c = £ and r c = 0. The results 
displayed in Fig. 3 show again that the order-parameter 
profile does not change qualitatively the core states. The 
overall shape of the spectra is determined by the strength 
of the supercurrent, which is proportional to the wind- 
ing number. For odd v, there is a branch crossing the 
Fermi energy at p = 0, while for even v there is no such 
branch. In addition, other branches cross zero energy at 
higher angular momenta. The profile of A(r) affects the 
core states with \p\ < \v\g. Except for small quantitative 
differences, the LDOS computed with and without the 
normal core arc therefore very similar, as for the singly- 
quantized vortex. 

From Figs. 1 and 2 one sees that the presence of the 
supercurrent is a necessary and sufficient condition to 
have low-energy states (while the normal core is not), 
and from Fig. 3 one sees that the topological frustration, 
which is only present for odd u, is a necessary condition 
to have low-energy and \ow-\p\ states, i.e. a zero-bias 
peak in the LDOS at the core center. In order to build 
a consistent explanation of the origin of the vortex-core 
states, it is thus necessary to understand the effect of the 
circulating superfluid on the BdG excitations. 



B. Supercurrent and pairing energy 

As already mentioned, the velocity field due to the su- 
percurrent changes the angular momenta of the electron 
and hole parts of the BdG excitations in opposite ways. 
As a result, a phase difference of is induced between 
the radial wave functions of the electron and hole. This 
phase difference is at the origin of the strong dependence 
of the vortex-core spectra on v. One possible way to solve 
the BdG equations is to treat the modulus of the order 
parameter perturbatively, while taking full account of the 
phase. For a vanishing modulus Eqs. (5a) decouple and 
the radial wave functions assume the simple form 

iP° ± (p)=A±J li ±»(gVT±^p), 

where J is the Bessel function. The eigenvalues e must be 
determined from the boundary condition at the border 
of a normalization disk of radius R, and they form a 
continuum for each value of p as R —* oo. The phase 
shift can be easily seen from the asymptotic behavior of 

^m(^); 

J m (x) ~ \f~ cos [ x ~ ( m + h) f] (a;»m), 

remembering that fi = n + ^ with integer n. The first- 
order perturbation theory in S gives the energies of the 
core states as 

e,=e + 2 f d P p5<j>)1> + { P )tf_{ P ). 
Jo 

Because of the phase shift, the pairing matrix element is 
qualitatively different for even and odd values of v. For 
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FIG. 4: Electron (ip+) and hole (ip-) BdG amplitudes for the 
lowest angular momentum in the singly- (top) and doubly- 
bottom) quantized vortex with r c — £. The integrand of 

the pairing matrix element, p tanh(p/-y/|f|)' I '''!/>+(p)' ! /'-(p); is 
displayed in the bottom panels. 



even v and small s, ip+ and tp?_ are either in phase or out 
of phase by n, and the integrand is thus either positive 
or negative, leading to a maximal pairing energy of order 
Aqo. For odd u, on the contrary, the integrand oscillates 
about zero and the resulting pairing energy is minimal. 
This mechanism is illustrated in Fig. 4, using the exact 
eigenstates rather than ip±. The solutions ip± as well as 
the integrand of the pairing matrix element arc displayed 
for the lowest allowed value of p and for winding numbers 
v = — 1 and —2. The numerical values of the kinetic and 
pairing energies, defined as 



£kin = (lp+\L + \lp+) 
£pair = 2{lp + \5\lp-), 



(ip-\L-\i>-) 



are also indicated. 

The analysis of the numerical results at all angular mo- 
menta shows that the energy eigenvalues in Figs. 1 and 3 
are dominated by the pairing energy. Due to the cancel- 
lation of the electron and hole contributions, the kinetic 
term remains small and has a weak p dependence. The 
structure of the eigenvalue spectra, and in particular the 
occurrence of low-energy states at high angular momen- 
tum for \u\ > 1, can thus be qualitatively understood 
by considering the evolution of the pairing matrix ele- 
ment with p. With increasing p, an additional phase 
shift appears between the radial electron and hole wave 
functions at short distances (as can be seen from the func- 
tions ip±(p), the phase shift always tends to asymp- 
totically, but it is a function of p for p < p/g). This 
has the effect of increasing the matrix element towards 
its maximum value of Aoo. For the lower branch in the 



v = — 2 spectrum in Fig. 3, there is a p > such that the 
total phase shift is close to — ^, and the matrix clement 
thus nearly vanishes. 

In the Appendix I derive an approximate solution for 
the bound states of the s-wave vortex. To lowest order 
in pv/gp c the eigenvalues are found to be 



-i^ + (2m+ ! + ,/)§ 
1 + gSp c 



(8) 



where to is any integer such that < A^ and p c is 
some effective core radius in units of £. For v = — 1, 
gS = -, and p c = 1, one recovers the famous Caroli-de 
Gennes-Matricon finding, 6 namely that for each p there 
is a unique solution within the gap (to = 0) with energy 

The term in Eq. (8) can be traced back to the phase 
difference discussed above (see the Appendix). Only for 
odd values of v does the term in parentheses disappear for 
some to, hence a branch of bound states crosses the Fermi 
level at p = 0, contributing to the low-energy LDOS in 
the core. In this case the gap between the lowest-energy 

A 2 2 

excitations in the core is A£^ =±1/2 = 2 p c (lr+2 Pc ) ' 
where p c itself is a function of v. The ^-dependence 
of AE^ = ±i/2 may be estimated by demanding that the 
"hole" in A(r) in Eq. (6) and in the model calculation of 
the Appendix have identical volumes, namely irpl. We 
then obtain p 2 c = 2\v\ dxx(l - tanh M x) w \v\ 2a with 
a = 0.78. Hence A£^ = ±i/ 2 is a decreasing function of 
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(y odd). (9) 



For even values of v there is always a large gap at p = 0, 
which is given by A£^ = o = A 00 7r/(1 + gSp c ) in this 
model. It should be noted that the presence (absence) 
of this large gap for even (odd) v is due to the ab- 
sence (presence) of the flux-induced topological frustra- 
tion, and does not depend on the order-parameter profile; 
however, the width of this gap does depend on the order- 
parameter profile through the effective core radius p c . 
Using p c w \v\ a and assuming the BCS relation gS = ^ 
holds, we find that A£^ =0 is also a decreasing function 

ofM, 



AE^o m A c 
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(10) 



in good agreement with the numerical results in Fig. 3. 



IV. THE ISOLATED D-WAVE VORTEX 

Within the BdG theory, there is an important qualita- 
tive difference between the core spectra of isolated vor- 
tices in ci-wave and s-wave superconductors. The energy 



spectrum of s-wave vortices is discrete (but looks contin- 
uous in experiments because of small interlevel spacing 
and thermal broadening), while numerical calculations 
suggest that the energy spectrum is continuous in the 
d-wave case, 19 although there are conflicting opinions in 
this respect. 20 A complete analytical solution would be 
useful to address this problem, but to my knowledge none 
has been published so far. The main difficulty comes 
from the nonlocality of the (i-wave gap operator: gauge- 
invariant generalizations of the lattice (i-wave gap to the 
continuum model are rather complicated. 21 In this sec- 
tion, the dependence of the core energy spectrum on the 
gap profile and vortex winding number is investigated nu- 
merically in the (i-wave case. The issue of bound versus 
extended core states is not directly relevant here. 

Instead of a continuum model I consider a two- 
dimensional nearest-neighbor lattice model. The vortex 
order parameter is taken to be 



±A(|i^|)cos(2T y )e"^ if \n -rj\=a 



otherwise, 



where rj 



denotes a lattice site, a is the lattice parame- 



r, = 



\Rij\(cos0ij, sin( 



0. 



a(cosTij, sin-Tjj), and A(r) is given by Eq. (6). The 
LDOS in the core is calculated using the Green's function 
formalism: 



■ImGuiE + iT). 



(11) 



The Green's function is given by the Gorkov equations 

G,» = GOH + ^G° fe HS H HG y H (12a) 
kl 

1 x - e ik-{ri-rj) 

= j^E^T (12b) 

k 

£;» = ^(r^GK-u^irurj) (12c) 

kl 

with iV 2 the number of k points and the dispersion. 
The details of Ek are not important in the context of this 
study. However, the presence of a van-Hove singular- 
ity in the gap region provides additional spectral weight, 
which might be unequally distributed among the various 
peaks in the spectra, thus complicating their interpre- 
tation. In order to avoid such difficulties, the simple 
nearest-neighbor form e fc = — 2t(cosk x a + cos k y a) — /i 
is used, with the chemical potential set to /i = t so that 
the van-Hove singularity does not influence the spectra 
in the gap region. Eq. (12) is solved by first computing 
Gij(u)) for a large system (N = 1024), taking advantage 
of the translational invariance in Eq (12b). Then the 
inhomogeneous terms Ejj(w) and Gij(u>) are calculated 
on a smaller M x M system (M = 51) with the vortex 
at the center. With this method the finite-size effects do 
not contaminate the free propagator, and a good spectral 
resolution can be achieved. The hopping integral is set 
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Eq. (6), r c = 2a 
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FIG. 5: Local density of states in the core of d-wave vortices 
with increasing winding numbers. The solid lines correspond 
to vortices having a normal core given by Eq. (6), and the 
dashed lines correspond to vortices without normal core. The 
case v = corresponds to a normal core without supercurrent. 
The arrows show the displacement of the peak maxima when 
the system size increases from M = 41 to M = 51. The 
bulk DOS is shown in grey, and provides a common scale 
to compare the spectra at different v. Note that the bulk 
coherence peaks are not located at E = ±Aoo due to the 
large value of the chemical potential fi = —t. 



to t = 5Aoo (a value typical for high-T c materials) and 
the energy broadening is T — Aoo/50. 

The LDOS calculated at the vortex center with r c = 2a 
and r c — and for various winding numbers is displayed 
in Fig. 5. The LDOS for a "normal core" of radius 2a 
without supercurrent is also shown and denoted as v = 0. 
One can see several striking similarities with the s-wave 
vortex. For \v\ ^ 1, the suppression of the order parame- 
ter has a very small effect on the LDOS in the core. This 
is confirmed by the case v = 0, which shows that a local 
suppression of A(r) is unable to induce low-energy states. 
(The small flat DOS at E = for v = and v = -2 is 
probably a finite-size effect. 22 ) For v = 0, there is a trans- 
fer of spectral weight from high to low energy at E ~ Aoo, 
resulting in two sharp states near the gap edges. Increas- 
ing the system size from M — 41 to M = 51, these states 
sharpen while moving to slightly lower binding energy, 
as indicated by the small arrows in the figure. At con- 
vergence, these states would presumably lie within the 
bulk gap, as in the s-wave case. For \u\ ^ 1, one observes 
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FIG. 6: LDOS along the nodal direction for d-wave vortices of 
increasing winding numbers. The curves show the LDOS at 
all sites between (0, 0) and (20, 20), and are offset vertically 
for clarity. 



a strong even/odd effect. Vortices with odd v have a 
zero-bias conductance peak (ZBCP), which is absent in 
even vortices. The ZBCP sharpens with increasing sys- 
tem size, but its energy is well converged at the system 
size considered. For \v\ ^ 2, there are additional states 
at high energy, which are the analogs of the high-energy 
states near \i = in the s-wave case (see Fig. 3). One 
can observe that the sensitivity of the peaks to system 
size increases with energy, suggesting that the states at 
higher energy are less localized. 

In the neighborhood of the vortex axis, the LDOS has 
many characteristics in common with the LDOS of s- 
wave vortices. In Fig. 6, one can see that in the case 
v = — 1 the ZBCP at the core center splits into two 
symmetric peaks at larger distances, very much like the 
high angular momentum electron-hole excitations in the 
s-wave case (compare Fig. 2). These excitations have 
been the focus of a recent self-consistent calculation. 23 
In the doubly-quantized vortex, the two peaks near the 
gap edges vanish slowly as the distance from the core in- 
creases, while a ZBCP develops, which is maximum at a 
distance r = 5y/2a. The latter again decomposes into two 
symmetric excitations at larger distances. The resulting 
LDOS is thus very similar to that of a doubly-quantized 
s-wave vortex (see Fig. 4 of Ref. 12), and can be read- 
ily interpreted on the basis of Fig. 3. A similar analysis 
shows that all of the features of the s-wave vortex with 
v = — 3 (Fig. 4 of Ref. 12) can be distinguished in the 
LDOS of the v = -3 d-wave vortex (Fig. 6). The LDOS 
exhibits some anisotropy around the vortex, but the dif- 
ferences between the nodal and antinodal directions are 



small. In particular, all the features discussed in Fig. 6 
are also present in the antinodal direction. 



V. PERTURBATIONS OF THE 
ORDER-PARAMETER PHASE 

The results of the previous sections suggest that the 
core states have a topological origin. An implication of 
this is that the states must not be suppressed by changes 
in the phase of the order parameter, as long as these 
changes preserve the topological defect. Conversely, the 
core states should be strongly affected by the interac- 
tion with an antivortex, since the latter annihilates the 
effect of the topological defect on orbits larger than the 
vortex-antivortex distance. In this section, I investigate 
the effects of some perturbations of the phase field on the 
core states in a <i-wave superconductor, starting with the 
perturbation induced by a nearby antivortex. 

The order-parameter for a vortex-antivortex pair was 
taken as the normalized product of the order parame- 
ters associated with each constituent (i.e. two vortices 
with winding numbers 1 and —1 and r c = 2a). The 
vortex is located at the origin and the antivortex at po- 
sition b = (—6, —b). In Fig. 7 is shown the LDOS at the 
vortex center and along the diagonal, in the direction 
opposite to the vortex-antivortex direction. The phase 
field around the vortex is also displayed. For b = 20, 
the LDOS is similar to the LDOS of an isolated singly- 
quantized vortex (Fig. 6). However, although the phase 
field is barely modified by the antivortex in the region 
where the LDOS is calculated, the latter is considerably 
broadened: the height of the ZBCP is only 71% of the 
corresponding height in the isolated vortex, while the gap 
integrated spectral weight is conserved within 2%. Note 
also that the ZBCP is already split at site (1, 1), unlike 
in the isolated case where this splitting occurs first at site 
(2, 2). Reducing the vortex-antivortex distance (b = 10), 
the ZBCP at site (0, 0) also splits. At the same time, the 
energy separation between the two states away from the 
core center is increased. Finally, at b = 2 the ZBCP dis- 
appears completely and the LDOS resembles the LDOS 
of a v — "vortex" (see Figs. 5 and 2). 

The spectra in Fig. 7 show that the formation of the 
core states is a nonlocal process: the LDOS at a partic- 
ular site depends on the phase winding in a region much 
larger than the core radius r c . In fact the height of the 
ZBCP at (0, 0) reaches 95% of its value in the isolated 
vortex for vortex-antivortex separations as large as ~ 70a 
(b = 50). This value sets an upper bound for the size of 
the aforementioned region. On the other hand, a clear 
splitting of the ZBCP at site (0, 0) occurs for b = 19 
or less, which sets a lower bound to 27a. These num- 
bers can be compared with the spatial extension of the 
core states. In the s-wave vortex the latter is typically 
£ w 2E F /(k F A oc ) for a state at E = [sec Eq. (A3)]. 
Using E-p — bt (the position of the Fermi energy with 
respect to the bottom of the band in the calculations), 
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FIG. 7: LDOS along the nodal direction (indicated by a thick 
gray line in the top part) for a vortex-antivortex pair in a 
d-wave superconductor for various vortex-antivortex separa- 
tions. The antivortex (white dot, only visible at the shortest 
distance) is at position b with respect to the vortex (black 
dot). The order-parameter phase in the region of the vortex 
center (excluding the d x 2_ y 2 symmetry factor for clarity) is 
represented in the form of a small arrow on each bond of the 
square lattice. 



FIG. 8: LDOS along the nodal direction for a d-wave singly- 
quantized vortex subject to random phase disorder. W char- 
acterizes the strength of disorder. The phase field in the core 
region is shown in the top part. The spectra corresponding to 
the ordered case (W = 0) are shown in grey for comparison. 



We note that the integrated LDOS is conserved within 
3% at all sites and for all disorder strengths considered. 



k F — n/(2a), and t = SAoo one obtains I « 32a, which 
falls within the bounds deduced from the numerical simu- 
lations. Thus, although the typical localization radius of 
the bound states might be different in d-wave and s-wave 
superconductors, 19 one may conclude that the LDOS in 
the vortex core "feels" the presence of an antivortex (or 
another vortex in a vortex lattice) at any distance shorter 
than the extension of the bound states. 

The sensitivity of the LDOS to perturbations of the 
phase which preserve the topological defect was studied 
by randomly disordering the phase field, i.e. replacing 
the order parameter 'F of the isolated vortex by *f/e lx 
where x is a random Gaussian variable centered at x = 
with variance W (FWHM). The resulting LDOS is dis- 
played in Fig. 8. At W = tt/8 and W = tt/4, the 
spectra are practically unchanged with respect to the or- 
dered case. No new structure is induced in the spectra 
up to the strongest disorder considered. At W = it/2, 
though, the ZBCP is reduced (but not broadened) and 
the energy separation between the two electron-hole ex- 
citations away from the core is decreased. This might be 
attributed to the loss of rotational symmetry, which leads 
to a stronger interaction between angular momentum 
eigenstates, and to a redistribution of spectral weight. 



VI. DISCUSSION 

The numerical results reported in Figs. 2 and 5 imply 
that the suppression of the modulus A(r) of the order 
parameter in the core can be safely ignored when dis- 
cussing the mechanism of formation of the vortex states 
in s- and d-wave superconductors. The main reason is 
the small magnitude of the superconducting gap with re- 
spect to the Fermi energy: a local suppression of A(r) 
is a very weak perturbation for excitations at the Fermi 
surface, which does not provide a substantial gain of pair- 
ing energy for a localized Bogoliubov excitation (unless 
the localization radius is smaller than the core radius, 
which in turn costs a large kinetic energy). Making use 
of Eq. (8) with v = 0, one can see that a local suppres- 
sion of A(r) can induce low-energy states in the limit 
Aoo/E-p 3> (fcpTc) -1 , which is never attained in the BCS 
superconductors for which Aoo/Ep w |(fc F r c ) _1 . 

Thus the origin of the bound states must be searched 
in the phase x( r ) °f the order parameter. Far from the 
vortex, the slow variation of x( r ) is known to induce a 
small density of states in the gap by the Doppler-shift 
effect. 24 The Doppler-shift approximation is not valid in 
the core region because there the superfluid velocity is 
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too large, but also because this approximation neglects 
the topological defect associated with the phase winding. 
A comparison of the right panels in Figs. 7 and 8 un- 
ambiguously shows that the microscopic details of x( r ) 
are much less relevant to the formation of the vortex 
LDOS than the "macroscopic" phase winding. This is 
further evidenced by the qualitative differences between 
even-f and odd-f vortices. The mechanism described in 
Sect. Ill B provides a natural explanation to these re- 
sults, since the interference between the electron and the 
hole, which is responsible for the formation of the bound 
states, is a global property of the wave function and is 
not destroyed by local perturbations of x(r). The is- 
sue of self-consistency in the order parameter (modulus 
and phase) is one of the major obstacles to a complete 
analytical calculation of the vortex-core LDOS. Never- 
theless, the results of the present study suggest that the 
LDOS calculated from the simplest possible "trial" order 
parameter, i.e. a uniform modulus A(r, 9) = Aoo and a 
geometric phase \(r, 9) = v9, should exhibit all the char- 
acteristic features of the exact solution. (An interesting 
exception might be the case where subdominant order 
parameters of different symmetries are brought about by 
the self-consistency.) 

In the presence of several topological defects, the wind- 
ing v in Eq. (1) depends on the number of defects con- 
tained in the path of integration. Based on the analysis 
in Section V, I tentatively argue that the structure of the 
vortex core LDOS is determined primarily by the average 
vorticity in the region occupied by the core states. This 
statement was checked by considering a variety of vortex 
and antivortex configurations. For example the effect of 
a v = 2 antivortex on the LDOS in the core of a v = — 1 
vortex was studied. The antivortex was located at posi- 
tion b = (—6, —b) as in Fig. 7. For small (b/a < 2) and 
large (b/a > 20) vortex-antivortex separations the LDOS 
was found to be close to that of an isolated vortex. Near 
b = 5a the spectrum was similar to the topmost case in 
Fig. 5, corresponding to v = 0. Finally in the other re- 
gions the LDOS exhibited intermediate shapes. In this 
geometry the vorticity is —1 in the region r < b\f2 and 
+1 in the region r > b^/2, so that the average vorticity in 
a region of radius R is v = 1 — 4 (-^) if b\[2 < R and —1 
if 6\/2 > R- Thus \v\ s=a 1 for small and large b, whereas 
\D\ w for 6- R/2. 

In the normal state of the superconducting cuprates, 
the behavior of the high-frequency optical conductivity 25 
and the presence of a large Nernst signal 26 have been 
generally attributed to vortex excitations, in the form of 
unbound vortex-antivortex pairs. Recently, Lee 27 argued 
that such vortices have to be "cheap" , i.e. must be free 
of core states in order to have a formation energy com- 
parable to the thermal energy ~ k^,T c . Experimentally, 
there is indeed convincing evidence from STM 28 ~ 31 and 
NMR 32,33 measurements that the core states are sup- 
pressed at T < T c in the cuprates, pointing to a failure 
of the simple d-wave BCS theory in the superconducting 
state. Lee pointed out that a straightforward extension 



of the BCS theory which includes phase fluctuations is 
unable to produce "cheap" vortices in the normal state. 
This argument relies on the assumption that the average 
vortex core LDOS in a vortex-antivortex soup is similar 
to the LDOS of an isolated vortex. As Fig. 7 shows, how- 
ever, the vortex core states can be very efficiently sup- 
pressed in a "BCS" vortex anti-vortex pair, even when 
the vortex-antivortex separation is much larger than the 
core radius. It is likely that the same phenomenon also 
occurs in more sophisticated models of the high-T c super- 
conductors. A systematic investigation of the core energy 
as a function of the vortex density in a state of fluctu- 
ating vortex-antivortex pairs would thus be helpful to 
elucidate the nature of the normal state in the cuprates. 

In this study, the effect of the order-parameter symme- 
try on the properties of the vortex states was left aside, 
and the emphasis was put on the similarities of the LDOS 
in s- and d-wave vortices, which illustrates the key role of 
the vorticity in both cases. It is well known that a sub- 
dominant order parameter can induce qualitative changes 
in the LDOS of BCS vortices. 19 Furthermore, the delicate 
question of the spatial extension of the core states in the 
d-wave case (exponentially localized or extended) was not 
addressed. These issues are a good motivation to search 
for a realistic analytical solution of the BdG equations 
for vortices in superconductors with nonlocal pairing. 



VII. CONCLUSION 

In the BCS theory, the superconducting state is char- 
acterized by a complex order parameter, with a modu- 
lus A(r) related to the superfluid density, and a phase 
X(t") related to the supercurrent. Although the electronic 
states bound to vortices have been generally attributed 
to the suppression of the superfluid density in the core 
region, it was found that the bound states are completely 
formed when this suppression is overlooked. In order to 
explain the formation of the core states, a new mecha- 
nism was proposed, which relies upon the influence that 
the vortex supercurrent exerts on the Bogoliubov excita- 
tions. In short, the Bogoliubov excitations localize onto 
vortices because the topology of the latter — not the nor- 
mal core — gives to the electron and hole components the 
possibility to avoid one another and to minimize their 
pairing energy. The spectral properties of vortices car- 
rying more than one flux quantum, as well as of vortex- 
antivortex pairs, are consistent with this idea. The pro- 
posed mechanism of localization results from the partic- 
ular form of the coupling mediated by the supercurrent 
between the electron and hole components, and is there- 
fore unique to superconductors. 
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APPENDIX A: APPROXIMATE SOLUTION FOR 
THE S-WAVE VORTEX 

In this appendix an approximate solution of the BdG 
equations for an isolated vortex in a s-wave superconduc- 
tor is derived. This solution differs somewhat from those 
given in Refs. 6 and 7, and applies to vortices of arbitrary 
winding number v. The purpose here is to highlight the 
role of v on the vortex-core spectrum. The radial wave 
functions is written as 



(Al) 



withq = 5 (C+#, k = g(C-h)^, C = Ul+5 2 W-£ 2 )^ 



and 7 = j (fl++fl-)^ 
they are related by 



. The B± are complex numbers and 



B- 
~B~4 




5 — ij 



*[(s+-s-)pc-f § ] 



(A4) 



Introducing the solutions f± and f± into Eq. (Al) and 
matching the wave-function and its derivative at p = p c 
one obtains 



B s [kJ s + {q-g s )Y s ](Y s +iJ s ) 

A s {q-g s )(J*+Y?)-{Y B J> a -J,Yiy 1 > 



where g± = g\f\ ± e and H n (x) = J n (x) + iY n (x) with J 
and Y the Bessel functions of the first and second kind, 
respectively. Eq (5) becomes: 



1 

g 2 



d 2 



l 



m 



dp 2 \p +2 H ± )dp 



f±+S(p)^.f T =0. 



(A2) 

Hn±»(g±p) is abbreviated as H±, and H± = dH±/dp. 
At this point I make two simplifications, (i) I consider a 
square gap profile: 5(p) — S9(p — p c ). As shown in the 
main text, the vortex-core spectrum is only weakly influ- 
enced by the details of 5(p). (ii) I use the approximations 




pv 9+ + 9- 
Pc 2g+9- 



, z fi(g+-g-)p 



These expressions were obtained from the asymptotic 
form of the Bessel functions, 

V TTX 

by expanding to first order in p/g and vjg. The ap- 
proximation (ii) is therefore inadequate for p ± ^ > g. 
Furthermore, in the second line, the term proportional 
to 1/p was evaluated at p = p c . This is justified for the 
calculation of the eigenvalues, since the latter are deter- 
mined by matching the wave- functions at p = p c . Note 
also that this approximation will break down if p c <C 1. 



The factor e ±w ^ in the expression for / H± is the con- 
sequence of the phase shift discussed in Sect. IIIB. With 
these simplifications Eq. (A2) can be solved exactly. I 
restrict the analysis to subgap states with e < 8. At 
p < p c the solution is simply f±(p) = A± where A± are 
real constants. At p > p c , the solution can be written as 



B±e l(q - 



g±){p-Pc) p -k(p-p c ) 



(A3) 



with s = ±, J± = Jn±%{g±Pc), and Y± = Y ll ±^(g±p c ). 
Eqs. (A4) and (A5) provide the relation between the co- 
efficients A + and A_, as well as the eigenvalue equation: 



tan 



(.9+ -9-)Pc 



IT' 

V 2 



-E + T)yJ 'S 2 + 7 2 



erj + + 7 2 — s 2 



(A6) 



where 



V = 



■ Ct"f 



7 — aS 



a = tan arg 



B_/A. 
B+/A, 



In Eq. (A6) g±, 7, and 77 are all functions of e. In order 
to make Eq. (A6) more tractable, we notice that ( rts | 
and g± ss g since S, 7, and e are small numbers. Taking 
C = 2 an( i 9± = 9 m Eq. (A5) we obtain B s /A s = 1 
and therefore a = 0. Eq. (A6) can then be solved to first 
order in e: 



E 



tan_1 te)+( 2m + 1 + i/ )§ 



1 + 



\gpc 



(A7) 



+ gSp c 



gpc 



(2m+l + i/)f 



1 + gSp c 

where m is any integer such that \E\ < Aoo and the sec- 
ond line is valid to first order in pv/gp c . I have checked 
numerically that the eigenvalues resulting from Eqs. (A6) 
and (A7) are in excellent agreement at all p. Eq. (A7) 
correctly accounts for the qualitative difference between 
odd-i/ and even-j/ vortices. Due to the approximations 
made, however, the validity of Eqs. (A6) and (A7) is lim- 
ited to the region \pi>\ < g, and the zero-energy states at 
\p\ > g in Fig. 3 are not reproduced. 
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